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We examine dynamic heterogeneities in a model glass-forming fluid, a binary harmonic sphere 
mixture, above and below the mode-coupling temperature T c . We calculate the ensemble indepen- 
dent susceptibility X<±(To) an d the dynamic correlation length ^4(r a ) at the a-relaxation time r a . 
We also examine in detail the temperature dependence of r a and the diffusion coefficient D. For 
higher temperatures we find that the standard Stokes-Einstein relationship, D ~ r^ 1 , holds, but at 
lower temperatures a fractional Stokes-Einstein relationship, D ~ t~ ct with a = 0.65, emerges. By 
examining the relationships between r a , D, x^( T a), and (,4(r a ) we determine that the emergence of 
the fractional Stokes-Einstein relationship is accompanied by a dynamic crossover from r a ~ e k2 ^ 4 

at higher temperatures to r a ~ e lKi at lower temperatures. 



I. INTRODUCTION 

The fundamental understanding of the drastic slowing 
down of a super-cooled liquid's dynamics and the even- 
tual formation of a glass remains an open question after 
decades of research. One of the challenges is finding and 
understanding the universal features of glassy systems 
related to the dramatic increase of the relaxation time 
r a (or viscosity 77) as the liquid is cooled. This challenge 
is complicated by the observation that several dynamic 
regimes exist. The dynamic regimes are usually identified 
by examining the temperature dependence of the relax- 
ation time T a or the diffusion coefficient D as a function 
of temperature T. The validity of different theories, and 
the range at which these theories are applicable, are of- 
ten tested through the fitting of r a to the predictions of 
the theories. 

The mode-coupling theory of Gotze and collaborators 
PQ in its original form predicted that upon supercooling 
the relaxation time diverges as a power law at the so- 
called mode-coupling temperature T c , t q ~ (T — T c )~ 7 . 
While it was quickly recognized that a true power-law 
divergence is absent, for many liquids there is a tem- 
perature range where a power law fits the temperature 
dependence of the relaxation time very well over about 
three decades of the change of the relaxation time [2]. 
For lower temperatures and longer relaxation times devi- 
ations from the power-law divergence are found in both 
experiments [3J and simulations (3H3- These deviations 
are usually interpreted as the signature of a new relax- 
ation mechanism involving activated dynamics [8]. 

Several interesting observations have been made about 
liquids' dynamics within the temperature range where 
mode-coupling like power law fits work well. This tem- 
perature range, which we will call hereafter the mode- 
coupling regime, is well accessible to simulational studies 
[llj . We note that many phenomena occurring in this 
temperature range have been interpreted using the so- 
called potential energy landscape [12]. Thus, this tem- 
perature region is also referred to as the landscape influ- 
enced regime [TTJ[T3]. At present, the significance of the 



fact that the simulational results can be interpreted or 
rationalized in terms of two rather different theoretical 
scenarios is unclear. 

One observation is the emergence of dynamic hetero- 
geneities upon approaching the mode-coupling temper- 
ature, and this feature has been extensively studied in 
simulations [3J rHH2"5] . Briefly, there are subsets of par- 
ticles that move either much faster or much slower than 
would be expected on the basis of a Gaussian distribu- 
tion of particles' displacements. Moreover, these sub- 
sets form clusters and the average size of these cluster 
increases upon supercooling. Extensive work has been 
done to determine quantitatively the characteristic size 
of dynamic heterogeneities. One way to approach this 
problem is to analyze the so-called four-point structure 
factor [T71 HHJ [53] , which is the structure factor calculated 
using the initial positions of the particles that have moved 
less than a certain distance during a specified time period 
(as discussed later, the distance is typically a fraction of 
the particle diameter and the time is often the relaxation 
time). An analysis of this structure factor results in a dy- 
namic correlation length £4, which we will refer to as the 
four-point dynamic correlation length. The length £4 is 
found to increase with decreasing temperature. Thus, the 
increase of the dynamic correlation length correlates with 
the increase in the relaxation time. It was found J3TJ [53] , 
that large systems need to be simulated for many relax- 
ation times for an accurate determination of £4. The large 
computational effort required to perform such studies re- 
sulted in that the size of dynamic heterogeneities was 
mostly examined above the mode-coupling temperature. 

Another interesting observation is the violation of the 
Stokes-Einstein relation [TS] H7[ [SB] . For a liquid in its 
stable thermodynamic state it is generally found that the 
Stokes-Einstein relation holds and thus the diffusion co- 
efficient is inversely proportional to the relaxation time, 
D ~ r^ 1 . For strongly supercooled liquids this relation- 
ship breaks down and is often replaced by the so-called 
fractional Stokes-Einstein relation, D ~ T~ a with a < 1. 
The Stokes-Einstein relation violation is often rational- 
ized in terms of the emergence of dynamic heterogeneities 
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Recent simulations also suggest that there is some sort 
of dynamical crossover between the onset of supercool- 
ing and the mode-coupling temperature. Bcrthicr et al. 
[2"5] reported interesting changes in the dependence of the 
relaxation time on system size for a variety of systems, 
including a binary harmonic-sphere mixture, at temper- 
atures around and above T c . They argued that these 
changes can be rationalized within the theoretical frame- 
work of the Random First Order Theory (which is usually 
assumed to include the mode-coupling theory) [30]. The 
latter approach allows for the competition between mode- 
coupling and activated dynamics. Berthier et al. argued 
that these dynamics are impacted by the finite system 
size in different ways and this results in the change in 
the finite size effects in the vicinity of T c . 

The results of the finite size effects study were fore- 
shadowed by an earlier dynamic correlation length study 
of the same binary harmonic-sphere mixture. Kob et al. 
[32] froze a set of particles randomly selected out of an 
equilibrium configuration to form a stationary wall, and 
examined the dynamics of the particles that were still 
allowed to move as a function of the distance z perpen- 
dicular to the wall. By analyzing the z dependence of 
the relaxation time, they extracted a dynamic correla- 
tion length £ dyn . If we need to distinguish £ dyn from 
other dynamic lengths we will refer to it as the point- 
to-set dynamic correlation length. Kob et al. found that 
this length initially increased, but then decreased around 

Initially it was believed that the decrease of £ dyn indi- 
cated a decrease in the characteristic size of dynamic het- 
erogeneities. However, by using a direct method to calcu- 
late a dynamic length scale £4, which, as discussed above, 
is associated with the size of dynamic heterogeneities, we 
showed that dynamic heterogeneities continue to grow 
below T c [55] . We observed, however, that there is a 
change in the relationship between r a and £4 within the 
mode-coupling regime. Upon supercooling r Q ~ e fc2 ^ 4 , 

but a crossover to r Q ~ e lS * was found at a tempera- 
ture slightly above T c . 

In this paper we provide some details of calculations 
reported in Ref. [33] and we thoroughly examine the 
arguments for the change in the relationship between r a 
and £4. The latter crossover remained hidden in previous 
simulations of a hard-sphere liquid [2T] [55] due to the 
somewhat narrower range of r Q and £4 accessible in those 
simulations. However, the crossover becomes apparent 
after examining the relationship between the relaxation 
time and the characteristic size and strength of dynamic 
heterogeneities combined with the analysis of the Stokes- 
Einstein relation violation. 

The paper is organized as follows. In Sec. [IT] we briefly 
describe the simulations performed for this work. In 



Sec. Ill we examine single-particle dynamics and the vi- 
olation of the Stokes-Einstein relation. In particular, we 
characterize in some detail the temperature dependence 
of the the relaxation time and the diffusion coefficient 



and compare them with the predictions of the mode- 
coupling theory, which allows us to identify the mode- 
coupling regime. While this identification is relatively 
straightforward, we found that there is some freedom 
and therefore some ambiguity regarding the identifica- 
tion of the mode-coupling temperature. We also look at 
two other possible functional forms for the temperature 
dependence of the relaxation time and the diffusion coef- 
ficient for temperatures where the mode-coupling power 
law fits are no longer appropriate. In Sec. [TV] we study 
the four-point structure factor. We use it to extract the 
so-called dynamic susceptibility and the dynamic corre- 
lation length. We examine the relationship of the latter 
quantities to the average dynamics, quantified by the re- 
laxation time. We show that the emergence of a frac- 
tional Stokes-Einstein relation found for low tempera- 
tures is consistent with a crossover in the relationship 
between the relaxation time and the dynamic correlation 
length. We summarize the results in Sec. [V] We dis- 
cuss variants of the procedure used to identify the mode- 
coupling temperature and provide details of the calcula- 
tion of the dynamic susceptibility in two Appendices. 



II. SIMULATION DETAILS 

We simulated a 50:50 mixture of harmonic spheres. 
The temperature and density dependence of this system's 
relaxation time has been investigated before [3H [35] . 
This system was also used by Kob et al. to study the 
temperature dependence of the point-to-set dynamic and 
static correlation lengths [32]. The interaction potential 
is given by 



V nm {r) = - [ I 



if r > a nm 



if r < ov 



(1) 



For our binary mixture (J22 = 1-4cth and a±2 = l-2aii, 
where the diameters are chosen to inhibit crystallization. 
Our results are presented in reduced units where an is 
the unit of length, \J ma\ x je is the unit of time, and 
10~ 4 e is the unit of temperature. We studied the system 
at a density p = N/V — 0.675. We ran simulations 
consisting of N = Ni + N 2 = 10, 000 particles for 30 > 
T > 6, N = 40,000 particles for 6 > T > 5, and N 
1 ()0. 000 particles for T = 5. At some temperatures we 
run simulations for several different N to check for finite 
size effects. No finite size effects were observed in the 
simulations presented in this work. We note that finite 
size effects reported by Berthier et al. [23] were found 
in systems much smaller than those investigated in our 
work. 

At all temperatures, we performed simulations in the 
NVE ensemble with a time step of 0.02 using the velocity 
Verlet algorithm. For T = 5, the lowest temperature 
studied, we found significant energy drift after several 
hundred million time steps. Thus, at T = 5.5 and 5 we 
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ran NVT simulations using a Nose-Hoover thermostat 
with a time step of 0.06. The simulations were run using 
LAMMPS code [Ml EI], which was modified to include 
the harmonic sphere potential. 

For T > 5.5 and every system size we ran four pro- 
duction runs of at least 100r Q (the a relaxation time, 
T a , is defined in Sec. III). Note that this includes four 



NVT simulations and four NVE simulations at T = 5.5. 
For T = 5 we ran four NVE simulations for 10r Q with 
N = 100,000 and 40,000 particles, and four NVT sim- 
ulations for 100T a with N = 40, 000 particles, and three 
NVT simulations for 100r Q with N = 100, 000. 

At every temperature we ran over 90r a in an NVT 
ensemble to equilibrate the system. To ensure that the 
T = 5 simulations were equilibrated, two of the 40,000 
particle simulations were equilibrated for 200t q and these 
simulations produced statistically the same results as the 
ones equilibrated for 100t q . Note that for T = 5 simula- 
tions of 100t q with a time step of 0.06 requires 1.4 x 10 9 
time steps, thus, including equilibration, the T = 5 runs 
ranged from 2.8 x 10 9 to 4.2 x 10 9 time steps. Addi- 
tionally, to check for equilibration we examined two and 
four-point time-dependent correlation functions for ag- 
ing. We found that the four-point susceptibility xa calcu- 
lated directly from the NVE and NVT simulations (for 
which we will later use symbols xa\nve and xa\nvt) is 
sensitive to aging as well as energy drift. We found no 
evidence of aging in our production runs. 

Finally, to check if the Nose-Hoover thermostat intro- 
duced any artifacts, we compared the results of the NVE 
and the NVT simulations at T = 5.5 for runs of 100r a in 
both ensembles. No measurable differences were found. 
Additionally, at T = 5 we compared an NVE simulation 
that was ten times shorter than the NVT simulations. 
The results (e.g. relaxation times and four-point struc- 
ture factors) agreed to within error. We should point out 
that poor choice of NVT simulation parameters can lead 
to artifacts; see Appendix B for more details. 



III. SINGLE PARTICLE DYNAMICS AND THE 
MODE-COUPLING CROSSOVER 

In this section we examine several aspects of single 
particle dynamics in super-cooled liquids. First, we ex- 
amine the overlap function, and we use it to obtain the 
a- relaxation time r a . We compare the temperature de- 
pendence of the overlap function with that of the self- 
intermediate scattering function. We also examine the 
mean-square displacement (5r 2 (i)), and calculate the dif- 
fusion coefficient D from the long time behavior of the 
mean-square displacement. We examine mode-coupling- 
like power laws and identify the mode-coupling regime, 
i.e. the temperature range where these power laws are 
a reasonable description of the data. In this section we 
discuss fits of the a-relaxation time and the diffusion coef- 
ficient that assume the same mode-coupling temperature 
for both quantities. In Appendix A we briefly discuss fits 



that do not use this constraint. Finally, we examine the 
Stokes-Einstein relation and find a crossover from the 
Stokes-Einstcin relation to a fractional Stokes-Einstein 
relationship, which works well at low temperatures. The 
temperature at which the relationship between D and r a 
changes from D ~ t^ 1 to D ~ T~ a plays an important 
role in the discussion in Sec. IIV1 In this section we al- 
ways analyze quantities pertaining to the whole system 
(i.e. involving both small and large particles). In Ap- 
pendix A we briefly discuss quantities pertaining to the 
individual components. 

We begin by examining the the overlap function 



F {t) = N~ l (X>n(*) 



(2) 



where w n (t) — 0(|r„(t) — r n (0)| — a), r n (t) is the position 
of particle n at a time t, is Heaviside step function, 
and, as discussed above, the sum runs over all particles of 
the system. The microscopic overlap function ^2 n w n (t) 
selects particles that have moved less than a distance a 
over a time t. We fix a = 0.3 so that the decay time 
of F Q (t) is approximately equal to the decay time of the 
self-intermediate scattering function F s (q; t) , 



F s (q-t)=N- 



V^ e -iq.(r„(t)-r„(0)) 



(3) 



for q corresponding to the peak of the total static struc- 
ture factor, q = 6.1. In Eq. |3| the sum runs over all 
particles of the system. The overlap function and the 
self-intermediate scattering function encode similar in- 
formation about the dynamics of our system. For a large 
system, calculating the former function requires signifi- 
cantly less computational effort. In addition, as discussed 
in the next section, we quantify dynamic heterogeneities 
in terms of correlations of the microscopic overlap func- 
tion pertaining to individual particles. For these reasons, 
we define the relaxation time in terms of the overlap func- 
tion. Specifically, the a-relaxation time r Q is the time at 
which the overlap function decays to e _1 , F (r Q ) = er 1 . 

Shown in Fig. [I] is F (t) for T = 30, 25, 20, 15, 12, 10, 
9, 8, 7, 6, 5.5, and 5. The features found in F (t) are 
similar to what has been seen in F s (q; t) calculated using 
only large particles [3U|3S] and the total F s (q; t) showed 
below. At high temperatures the overlap function agrees 
with what is expected from a Gaussian distribution of dis- 
placements (note that the high temperature limit of the 
overlap function can be expressed in terms of the error 
function and is different from a simple exponential). At 
low temperatures a plateau develops followed by a slow, 
broad decay to zero. The low-temperature final decay 
is well described by a stretched exponential A e"W T °^° 
with 0.6 > j3 > 0.52 for temperatures between 8 and 5. 

The large number of particles used in our simula- 
tions makes the calculation of the intermediate corre- 
lation functions computationally expensive. Thus, we 
calculated the total self-intermediate scattering function 
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FIG. 1: The overlap function F (t) for T = 30, 25, 20, 15, 12, 
10, 9, 8, 7, 6 , 5.5, and 5, listed from left to right. The dashed 
line is a stretched exponential, Aoe'^^"^" , fit to F a (i) at 
T = 5 where /3„ = 0.53. 




FIG. 2: The self-intermediate scattering function F s (q;t) for 
q = 6.1 and T = 30, 25, 20, 15, 12, 10, 9, 8, 7, 6, 5.5, and 
5, listed from left to right. The dashed line is a stretched 
exponential, A s e~' t//Ts - ) 3 , fit to F 3 (q; t) at T = 5 where fi s = 
0.55. 



F s (q;t) defined in Eq. ([3| for a single wave- vector only. 
For F s (q;t) calculated using all particles, we used the 
wave- vector corresponding to the peak of the total struc- 
ture factor, q = 6.1. Shown in Fig. [2] is F s (q;t) for 
T = 30, 25, 20, 15, 12, 10, 9, 8, 7, 6, 5.5, and 5. At 
high temperatures the self-intermediate scattering func- 
tion exhibits an approximately exponential decay. As 
was the case with the overlap function, at lower temper- 
atures the self-intermediate scattering function develops 
an intermediate plateau followed by a slow, broad decay. 
The low-temperature final decay is well described by a 
stretched exponential A s e~ { - t / Ts ' ) ' is with 0.77 > j3 s > 0.55 
for temperatures between 12 and 5. 



the inverse temperature. 



In Fig. [3] we compare the temperature dependence of 
the a-relaxation time defined above and the more usual 
a-relaxation time T as defined in terms of the decay of the 
self-intermediate scattering function, F s (q;T as ) = e _1 . 
Both relaxation times exhibit very similar temperature 
dependence. In the same figure we also show characteris- 
tic times obtained from the stretched exponential fits to 
the final decay of the overlap function, t d and the self- 
intermediate scattering function, r s . These times can 
only be determined at low enough temperatures. Again, 
they exhibit a temperature dependence similar to that of 
the a-relaxation time. 

In Fig. [4] we show the temperature dependence of the 
amplitude of the stretched exponential fit to the self- 
intermediate scattering function A s and the stretching 
exponents /3 and /3 S . The former quantity, A s , is signifi- 
cant in that it may be used to identify the mode-coupling 
regime. Specifically, the mode-coupling theory predicts 
that the self-intermediate scattering function develops a 
plateau which extends to longer and longer times upon 
approaching the mode-coupling transition. The height of 
this plateau is predicted to be temperature independent. 
While extracting the plateau height from data showed in 
Fig. [2] is somewhat delicate, Weysser et al. [38] showed 
that the amplitude of the stretched exponential fit agrees 
with a value obtained from a more sophisticated proce- 
dure. 

Fig. |4ja) shows that the amplitude of the stretched 
exponential fit is approximately temperature indepen- 
dent down to T = 7 and then starts increasing. This 
would suggest that the mode-coupling regime ends at 
this temperature. While the data showed in Fig. Qa) 
are quite suggestive, we feel that this issue needs to be 
investigated further before making a strong claim that 
the mode-coupling regime ends at T = 7. Therefore, in 
the fits described below we use two different temperature 
ranges, 12 > T > 7 and 12 > T > 6. We note that the 
upper limit of the mode-coupling regime is determined by 
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FIG. 4: (a) The amplitude of the stretched exponential 
fit to the long-time decay of the self-intermediate scattering 
function as a function of the inverse temperature. Accord- 
ing to Ref. 38 , A s is a good estimate of the intermediate 
time plateau of F a (q;t). (b) The stretching exponents /3 and 
I3 S obtained from the stretched exponential fits to the long- 
time decay of the overlap functions and the self-intermediate 
scattering function, respectively, plotted versus the inverse 
temperature. 

the fact that above T = 12 a plateau in either the over- 
lap function or the self-intermediate scattering function 
cannot be determined. The latter observation agrees rea- 
sonably well with value of the onset temperature of slow 
dynamics which has been found to be around T = 13 by 
Berthier et ol. [22] . Finally, we note that Fig. |4jb) shows 
that at low temperatures the stretching exponents ob- 
tained from fitting the final decays of F (i) and F s (q;t) 
to stretched exponentials converge. We note that the 
overlap function can formally be obtained through an in- 
tegral of the self-intermediate scattering function over all 
wave- vectors, 

Fo(t) = J j0^f(q;a)F s (q;t) (4) 

where f(k; a) = Aita 2 ji(ka)/k with j\ denoting a spheri- 
cal Bessel function of the first kind. Eq. Q together with 
the low-temperature convergence of the stretching expo- 
nents suggests that the characteristic relaxation time t s 
obtained from the stretched exponential fit to the long- 
time decay of F s (q; t) and the stretching exponent (3 S are 
wave- vector independent, at least for the wave- vectors 
making the dominant contribution to the integral in Eq. 
Q. To investigate this further one needs to evaluate 
and analyze the self-intermediate scattering function for 
a range of wave- vectors. 

Another measure of single particle dynamics is the 
mean-square displacement 

(5r*(t))=N-UY / K(t)-r n (0)} 2 ), (5) 



FIG. 5: The mean square displacement (8r 2 (t)) for T = 30, 
25, 20, 15, 12, 10, 9 8, 7, 6, 5.5, and 5, listed from left to right. 



which is shown in Fig. [5] We clearly see a short time 
ballistic regime where (Tr 2 (t)") = 3Tt 2 , and a diffusive 
regime where (i5r 2 (t)) = 6Dt. Furthermore, at low tem- 
peratures at intermediate times a plateau region emerges 
between the ballistic and diffusive regimes. This plateau 
is associated with the cage effect where particles are 
trapped inside a cage of neighboring particles. How the 
particles escape this cage is one of the fundamental ques- 
tions of the dynamics of glass-forming liquids. 

We determine the diffusion coefficient D from the 
asymptotic, linear in time, dependence of the mean- 
square displacement, D — \im t ^oo (8r 2 (t)} / (6t) . In 
Fig. [6] we show the dependence of the a-relaxation time 
and the diffusion coefficient on the inverse temperature. 
Lines in this figure indicate fits discussed in the next two 
paragraphs. 

According to the mode-coupling theory there should 
be one temperature T c at which the ergodicity is broken 
and one critical exponent 7. To allow for better quality 
fits we followed Kob and Andersen [5] [B] and we fitted the 
a-relaxation time data and the diffusion coefficient data 
to mode-coupling-like power laws t q = a(T — T c ) -7 and 
D = b(T — Tc) 1 . In other words, we imposed one mode- 
coupling temperature T c common to both quantities but 
we allowed for different critical exponents. As mentioned 
above, we used two different temperature ranges. Fits 
using the 12 > T > 6 result in T c = 5.1 ± 0.1, Y = 
2.9±0.1 and 7° = 2.1±0.1. Fits using 12 > T > 7 result 
in T c = 5.6 ± 0.1, j T = 2.4 ± 0.1 and 7° = 1.9 ± 0.1. All 
these fits are shown in Fig. [6] as dashed lines. 

Also shown in Fig. [6] are two other commonly used 
fits. Solid lines show r Q = T Q e k ° T ~ 2 and D = D e~ d ° T ~ 2 
and dotted lines show r a — T v e klJ ( T ~ T °^ and D = 
£>„e- d " (T ~ T ° r . Both sets of fits were performed for 
T < 8. Note that the latter set of fits imposes the same 
Vogel-Fulcher temperature T for both the a-relaxation 
time and the diffusion coefficient. This set of fits results 
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FIG. 6: The relaxation time r a and the inverse diffusion coef- 
ficient Z? _1 versus 1/T. The dashed lines are mode-coupling 
power law fits for two different fit ranges: 12 > T > 6 and 
12 > T > 7. The solid lines are fits to r a = T e k ° T ~ 2 and D = 
D e~ d ° T and the dotted lines are fits to r a = T v e kv{T ~ To) 
and D = D„e" tMT " T ° r 1 where T a = 1.7. The latter fits were 
performed using T < 8 data. 

in T = 1.7. We shall emphasize that either e BT 2 or 
e A(T-T a ) g^. g (jggcribe very well the low temperature 
data. Thus, our simulational results cannot exclude ei- 
ther a T = or a finite temperature divergence of the 
a-relaxation time. 

We note that the mode-coupling temperature obtained 
in this work from fits using the range of temperatures 
12 > T > 6, T c = 5.1 ± 0.1, is slightly smaller but within 
error of the T c = 5.2 reported by Kob et al. [32 . The 
latter value was obtained from fitting the a-relaxation 
times obtained from two separate self-intermediate scat- 
tering functions involving small and large particles. In 
Appendix A we show that the fit of the a-relaxation ob- 
tained from the overlap function involving all particles 
results in a value comparable to that of Kob et al. In 
the same Appendix we also discuss more elaborate fit- 
ting procedures using relaxation times and diffusion co- 
efficients pertaining to small and large particles. The dis- 
cussion above and in Appendix A suggests that instead 
of the single mode-coupling temperature T c we should 
introduce a range of mode-coupling temperatures or a 
T c range. We propose to identify the temperature range 
5.6 — 5.1 as the T c range for the system studied in this 
work and we will indicate this range in the figures rather 
than showing a single mode-coupling temperature. 

The above analysis shows that our results are very sim- 
ilar to what is experimentally observed for super-cooled 
liquids [2] . There is a range of temperatures where mode- 
coupling like power laws provide a reasonable description 
of the data. While the mode-coupling regime is reason- 
ably well defined, small changes in the range of temper- 
atures used for fitting result in somewhat large changes 
of the mode-coupling temperature. Importantly, neither 
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FIG. 7: Stokes-Einstein decoupling. For T > 12, D ~ r" 1 , 
but for T < 8 we find that the fractional Stokes-Einstein 
relationship D ~ r^ ' 65 provides an accurate description of 
the data. The thin solid vertical line indicates the relaxation 
time corresponding to temperature T = 8. 



the relaxation time diverges nor the diffusion coefficient 
vanishes at the T c inferred from the fits, and any mode- 
coupling-like transition is avoided. Furthermore, the ex- 
ponents associated with the power-law fits for r Q and for 
D differ from each other, which is contrary to the pre- 
dictions of mode-coupling theory. 

The difference in the mode-coupling-like power laws' 
exponents is not surprising when one realizes that the 
Stokes-Einstein relation is violated, see Fig. [7j For 
high temperatures we find that Stokes-Einstein relation 
is obeyed, D ~ t^ 1 till about T = 12. In contrast, 
for the lowest temperatures we find that a fractional 
Stokes-Einstein relation works well, D ~ t^ 65 . This 
fractional Stokes-Einstein relation emerges for T < 8, 
which is inside both temperature ranges used for the 
mode-coupling-like power law fits. Thus, the exponent 
for the fit to t q and D would not expected to be the 
same. Also, we note that there is a temperature range, 
12 > T > 8, where the standard Stokes-Einstein relation 
is violated, and the fractional Stokes-Einstein behavior 
has yet to emerge. 

In summary, for temperatures lower than what is well 
described by the power laws, a new behavior emerges 
which is well described by an exponential divergence of 
the relaxation time and an exponential vanishing of the 
diffusion coefficient. The nature of this divergence is un- 
der debate since different functional forms describe many 
different glass formers well (see e.g. Refs. [2 [39] and ref- 
erences therein). In the following sections, we will ex- 
plore the relationship between the liquid's dynamics and 
dynamic heterogeneities. 
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IV. DYNAMIC SUSCEPTIBILITY AND 
CORRELATION LENGTH 

In the previous section we examined different aspects 
of the single particle dynamics within the model glass- 
forming liquid studied here. We observed a change in 
many aspects of the dynamics as the temperature de- 
creases. In particular, the final relaxation at low temper- 
atures is broadened and is well described by a stretched 
exponential. Stretched exponential relaxation can result 
from different populations of particles where some popu- 
lations relax much faster or much slower than would ex- 
pected on the basis of a Gaussian distribution. To look 
for these subsets of particles we examine the self- van Hove 
correlation function, G s (Sr;t) = (5r — [r n (t) — r n (0)]}. 
Since short time ballistic and long time diffusive motion 
both result in Gaussian distributions of G s (5r; t), the ex- 
istence of fast and slow sub-populations may be found 
by comparing G s (5r;t) to a Gaussian distribution with 
the mean-square displacement equal to that calculated 
from the simulation. Following a procedure suggested 
previously [401 141) , we present here the probability distri- 
butions of the logarithm of single-particle displacements 
which is easily obtained from the self-van Hove corre- 
lation function, P[\og w (Sr);t] = ln(10)47nrr 3 G s (c5r, t). 
The advantage of looking at this probability distribution 
is that if the van Hove function is Gaussian, the shape of 
P[log 10 (Sr);t] is independent of time. 

Shown in Figs. [8] and [9] are the calculated probabil- 
ities of the logarithm of single particle displacements 
at the a-relaxation time for T = 20 and T = 5 com- 
pared to what would be obtained from a Gaussian van 
Hove function with the same mean-square displacement 
(5r 2 (r a )), calculated using the smaller and larger par- 
ticles. The Gaussian distribution describes the results 
well for T = 20, before the onset of supercooling. In 
contrast, at T = 5, which is below the mode-coupling 
temperature, we observe several peaks, which suggests 
an activated, hopping-like relaxation mechanism. While 
Figs. [8|9] clearly indicate the existence of sub-populations 
of slow and fast particles, they cannot reveal spatial cor- 
relations of these sub-populations of particles. 

To examine the spatial correlations of slow particles we 
followed previous work [TTJ, [JHJ 0UJ H2] and examined the 
four-point structure factor 



K(0)-r m (0)] 



AO) 



(6) 



Si(q; t) is the structure factor calculated using the parti- 
cles that have moved less than a distance a during a time 
t. It measures correlations between microscopic overlap 
functions pertaining to individual particles. For t = 0, 
S/i(q;t — 0) is equal to the static structure factor. At 
later times there is an enhancement of the low q values 
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FIG. 8: Comparison of P[\og 10 (5r);r a ] for the smaller parti- 
cles at T — 20 and T = 5 (solid lines) with what would be 
expected from a Gaussian distribution (dashed lines) with the 
same mean-square displacement. 
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FIG. 9: Comparison of P[log 10 (<5r); r a ] for the larger particles 
at T — 20 and T — 5 (solid lines) with what would be ex- 
pected from a Gaussian distribution (dashed lines) with the 
same mean-square displacement. 



which indicates that the slow particles, i.e. the particles 
that have moved less than a over the time t, form clus- 
ters (see Fig. [l0]for the four-point structure factor at the 
a-relaxation time). To examine these clusters we study 
the dynamic susceptibility Xi(t) = lim g — s-o S^iq; t), which 
is related to the number of dynamically correlated parti- 
cles, and the dynamic correlation length ^i(t), which is a 
measure of the size of the slow regions. 

To aid in fits which we use to determine the correlation 
length the we calculated 



Xi (t) = N- 



^2 w n(t)w m (t) ) - ( y^ y w n (t) 



(7) 




100 




FIG. 10: The four-point structure factor at some represen- 
tative temperatures. The solid lines are fits to the Ornstein- 
Zernicke equation X4(t q )/[1 + (£(t q )</) 2 ] for q < 1.5/^4(r a ). 



from simulations and then used it at t = r a as a fit point 
at S±{q = 0; r a ). We note that Xi(t) as given by Eq. Q 
is equal to the small wave-vector limit of the four-point 
structure factor lim 9 _>.o Si(q; t) only if the average in Eq. 
([7]) is performed in the ensemble that does not suppress 
global fluctuations jHHS]. Thus, one can not directly 
calculate it in the NVE and NVT simulations since en- 
ergy and/or particle number fluctuations are suppressed 
in these simulational ensembles. However, Xa{^) can be 
calculated by taking into account these suppressed fluc- 
tuations as described in Refs. [2IJ [22l 02] ■ Details of this 
method for the system studied here are given in Appendix 
B. For simplicity and to aid in the following discussion we 
introduce the following notation for a contribution due to 
energy fluctuations Xi{t)\sT = kBT 2 x^(t) / c v where ks is 
Boltzmann's constant, c„ is the constant volume specific 
heat per particle, and xr(t) = dF {t)/dT. In addition, 
we introduce the following notation for a contribution 
due to particle numbers fluctuations Xi(t)\sN- The ex- 
pression for the latter quantity is rather long and it is 
given in Appendix B. As discussed before [2T1 1221 |4"2H4"4] 



Xi(t) = Xi(t)\NVE + Xt(t)\sT + Xi{t)\6N (8) 

where Xi(t)\NVE is X±(t) given by Eq. Q calculated 
in an NVE simulational ensemble. Likewise we denote 
Xi{t)\NVT for Xi{t) calculated in the NVT simulational 
ensemble. At T = 5.5 we checked that 



Xi(T a )\NVT = Xi{ T a)\NVE + Xi(T a )\sT 



(9) 



by calculating all terms in Eq. ([9| separately. For 
T = 5, due to energy drift in the NVE simulations, 
caused by the very long simulation time (over 4 bil- 
lion time steps) required, we used Xi( T a)\NVT for the 
sum Xi{ T a)\NVE + X4:( T a)\sT- It is important to note 
that poor choices of simulation parameters can result in 
Xi{ T a)\NVT j= Xi{ T a)\NVE + Xii T a)\sT due to too strong 



FIG. II: The different terms which represent the main con- 
tributions to the ensemble independent susceptibility Xi( r a)- 
The circles are X4( t o)\nve and are directly calculated in an 
NVE simulation (open circles) except for at T = 5 (filled 
circle) where we used Xi{ T a)\NVE = X4(T a )| NVT — X^{ t ci)\st- 
The dashed line through the circles serves to guide the eyes. 



or too weak of a coupling to the temperature bath; see 
Appendix B for more details. 
Shown 



Fig. II 



are Xii T a)\NVE (open circles), 
X4( T a)\sr (squares) and Xi(T a )\sN (triangles) for T < 15. 
At high temperatures, Xa{ t o)\nve is the largest term. 
For T < 10, Xi(T a )\NVE grows slower than Xi(T a )\sN 
and for T < 6 the latter term becomes the largest con- 
tribution to Xi( T a)- For T < 8 both Xa{ t o)\sn and 
X4( t q)I<5t grow as T~ 4 with decreasing temperature. We 
note that T ~ 8 is the temperature at which the frac- 
tional Stokes-Einstein behavior emerges for this system, 
Fig. [7] At T — 5 we could not perform a long enough 
NVE simulation without significant energy drift, and 
a calculation Xi( T a)\NVE for short simulations results 
in a very large uncertainty. However, by looking at 
Xi( T a)\NVT ~ Xi{ T a)\sT it appears that Xi{ T a)\NVE at 
T = 5 decreases. If the above discussed trends continue 
to lower temperatures, then Xi( T a) should also grow as 
T~ 4 with decreasing temperature. 

Having 64(9; T a ) for q > from direct simulations 
and using Xi( T a) from the above described procedure as 
Si(q = 0;t q ), we now fit the four-point structure factor 
to an Ornstein-Zernicke function 

Xi(T a ) 



1 + (Ur a )q) 2 



(10) 



We determined previously [211 122] that these fits are ro- 
bust if the range of wave- vectors used in the fits is re- 
stricted to q < 1.5/^4(r Q ). The Ornstein-Zernicke fits 
are showed as solid lines in Fig. [TUj We note that the 
fits result in values of Xi{ T a) that are slightly different 
from those determined by the procedure described above. 
These differences are seen in Fig. [10] as differences be- 
tween data points at q = and the positions of the solid 
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FIG. 12: Comparison of different dynamic correlation 
lengths as a function of 1/T. The circles are £,4,(r a ) calcu- 
lated in this work, the left triangles are £c yn and the squares 
are £s yn . The latter two lengths are the dynamic point-to-set 
lengths described in Ref. [32]. The solid vertical line shows 
the maximum of £ dyn and the dashed lines gives the range of 
mode-coupling temperatures that are consistent with mode- 
coupling like power law fits to r a and D. 



lines at q = 0. In the remainder of the paper, we will 
plot and discuss Xi{ T a) as determined from the Ornstein- 
Zcrnicke fits. 



In Fig. 12 we compare the temperature dependence of 
the four-point dynamic correlation length determined in 
this work and the point-to-set dynamic lengths of Kob 
et al. These lengths have a similar magnitude at high 
temperatures. The four-point length £4 increases with de- 
creasing temperature faster than either of the two lengths 
determined by Kob et al. Most importantly, £4 increases 
with decreasing temperature for all temperatures exam- 
ined in this work whereas the lengths determined by Kob 
et al. peak around T = 6 and then decrease with de- 
creasing temperature. We note that the latter behavior 
correlates with the non-monotonic temperature depen- 
dence of the finite size effects investigated by Berthier 
et al. |29j . Curiously, our four-point dynamic correla- 
tion lengths are significantly longer than the point-to-set 
static correlation lengths determined by Kob et al. This 
contrast between four-point and static point-to-set corre- 
lation lengths is also present in other systems (compare, 
e.g. dynamic lengths determined in Ref. [5U] and static 
lengths of Ref. g5]). 

Next we study the relationship between Xi( T a) and 
£4 (fa) j which is shown in Fig. |13[ As in previous work 
[21] [22] . we find that Xi{ T a) ~ £ 4 ( r a) 3 when £;4(t q ) is 
greater than approximately 3, which occurs at a temper- 
ature of approximately 9. This result suggests compact 
domains of slow particles. We do, however, caution on 
obtaining a correlation length from ~ [xiit)] 1 ^. 

It was observed in previous work, Ref. that 
reaches a plateau and remains constant as a function of 



FIG. 13: The relationship between X4,( T a) and £4(7"^) where 
Xi( T a) and (,4(tc) are both obtained from the Ornstein- 
Zernicke fits to S4(q; r a ). The straight line is a fit to Xi( T a) ~ 
£ 4 (T a ) 3 , and describes the data well for T < 9. The inset is 
^(to,) 3 / 2 versus T~ 2 , which describes the data well over the 
whole temperature range. 



time during a period when Xi(^) decreases. Thus, to 
compare £4^) at different temperatures, different times, 
under confinement, etc. one should not use the relation- 

ship£ 4 W~[x 4 W] 1/3 . 

We now examine the relationship between r a and 
£&(T a )- Previously we observed for a hard-sphere glass- 
forming system that r Q ~ provided a good descrip- 
tion of the results over the full range of volume fractions 
studied in Refs. [2TTI221. 



In Fig. 14 we show r Q versus £ 4 (r a ) for the harmonic 
spheres. For T > 8 we find that r Q = r 2 e fc2 ? 4 , but for T < 
8 we find a faster growth of the a-relaxation time with the 

, >.3/2 

dynamic correlation length, r a = T\e ^ . We note that 
the Random First Order theory j31j predicts that r a ~ 
e^ 7 where £ s is a time-independent static correlation 
length characterizing the size of correlated regions, and 
this functional form is consistent with our results. 

Exponential correlations shown in Fig. [14] combined 
with power law relations between the diffusion coefficient 
and the a-relaxation time, D ~ t^ 1 and D ~ r^ 0,65 
(see Fig. [7|, imply similar exponential correlations be- 
tween the diffusion coefficient and the dynamic correla- 
tion length. In Fig. [l4]we also show 1/D versus £ 4 (tq,) 
for the harmonic spheres. We find that D = I^e 2 * 4 
fits the data well for all temperatures except T = 7. 
However, the T < 8 data are also compatible with 



; 3/2 



D = D\e~ aiti ^ . Simulations for somewhat lower tem- 
peratures, which do not seem feasible at present, are 
needed to shed some additional light at the relation be- 
tween the diffusion coefficient and the dynamic correla- 
tion length. 

We note the internal consistency of our low temper- 
ature correlations. If the observed trends continue then 
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FIG. 14: The a relaxation time r a (circles, left axis) and 
D^ 1 (squares, right axis) as a function of ^4,(T a ). The solid 

lines are fits to r„ = T\e klti ' and D = D\e~ dlii , while the 
dashed lines are fits to r a = r 2 e k2U and D = D 2 e~ d2ii . 



for low temperatures X4,(r a ) ~ Xi( T a)\&T +Xi( T a)\sN and 
Xi( T a) ~ T~ 4 . The last dependence, taken together with 
the correlation X4( t q) ~ implies that at low temper- 
atures £ 3/2 - T' 2 . We find that £ 3/2 - T" 2 describes 
the data well for the whole temperature range studied in 
this work even for temperatures wher e X A( T a) does not 
increase as T~ 

,k T~ 



-4. 

-i-2 



see the inset to Fig. 13 



Since we find 

that r a ~ e fc ° J provides a good description of the data 
for T < 9 (see Fig. ^J, one would expect that t q ~ e fcl ^ 3/2 
to also provide a good description for these temperatures. 

Our data is also consistent with a divergence of £4 ~ 
(T — To) -2 ' 3 , but to obtain a To consistent with the fits 
to T a and D we must introduce an offset and fit e.g. 
£4 = (o$ + bt/(T - T )) 2 / 3 . A similar sort of offset is also 
needed for fits to Xi{ T a)- 

In summary, in this section we examined the relation- 
ships between £4(7^), r a and D. Our results are consis- 
tent with r a = r 2 e k2U and D = D 2 e~ d2U for T > 8 and 



then a crossover to r a = T\e kl ^ 1 and D = Die~ dl ^' for 
T < 8. We also found that this relationship is consistent 
with the temperature dependence of the correction terms 
to Xi{ T a) and the observed fractional Stokes- Einstein re- 
lation examined in the last section. 



c 3/2 



V. CONCLUSIONS 

We examined dynamic heterogeneities above and be- 
low the mode-coupling temperature and found evidence 
of a dynamic crossover associated with the dynamic het- 
erogeneities. To identify the mode-coupling temperature 
we fit the a-relaxation time and the diffusion coefficient 
to mode-coupling like power laws. We obtained a slightly 
different mode-coupling temperatures T c depending on 
the range of temperatures used for fitting. We found that 



mode-coupling like power law fits to r Q and and D give 
different exponents for this system. This is consistent 
with the onset of the Stokes-Einstein relation violation 
inside the mode-coupling regime. Specifically, we found 
that D ~ t~ ' 65 beginning around T = 8. However, an 
exponential divergence of r a and D with the same form 
for both r a and D is consistent with a fractional Stokes- 
Einstein relation, and we find that two commonly used 
forms e BT and e ■M T ~ T °) are both consistent with our 
data. 

We examined and characterized dynamic hetero- 
geneities from T = 20 down to T = 5, which includes 
temperatures above the onset of slow dynamics, T « 13 
[29], and below the mode-coupling temperature range, 
T c m 5.6 — 5.1. We calculated the four-point suscepti- 
bility Xi{ T a) and the correlation length ^4(t q ) by fitting 
the four-point structure factor Si{q;r a ) to an Ornstein- 
Zernicke function. Through the relationships between 
D, r a , and £4 we found a dynamic crossover that had 
remained hidden in a previous hard sphere simulation 
pTl 122] due to a smaller range of relaxation times. We 
found that r Q ~ e k ^U f or T > 8, which includes the first 
decade of slowing down after the onset of slow dynamics 
at T w 13, then t q - e kl &' 2 for T < 8, where T = 8 
is the same temperature where we find the onset of the 
fractional Stokes-Einstein relationship. In contrast, we 
found that D ~ e ~ d2 ^ 1 fits our data well over the whole 
temperature range. We note that to be consistent with 
the fractional Stokes-Einstein relationship one would cx- 

pect D ~ e 144 for T < 8 and we found that this is 
also consistent with our data. Simulations at lower tem- 
peratures are needed to examine the relationship between 
D and £. 

For T < 8 there is some interesting behavior of dif- 
ferent contributions to the four-point dynamic suscep- 
tibility X4:( T a)- In our simulations we directly calculate 
the constant volume and energy contribution Xa{ t o)\nve 
or the constant volume and temperature contribution 
Xi( T a)\NVT, and to obtain the lim^o of S4(q;r a ) we 
calculate fluctuations suppressed in the NVE or NVT 
simulations. These correction terms are related to the 
derivatives of the two point correlation function Xx (t) = 
dF (t)/dx. We find that Xp and Xt both grow as T~ 4 
starting around T = 8; again around the same temper- 
ature where the fractional Stokes-Einstein relationship 
emerges. 

We see evidence that Xi{ T a)\NVE levels off around 
T = 6 and decreases around T = 5. The similarity of the 
temperature dependence of X4( T a)l-/vv.E and the point- 
to-set dynamic correlation length £ dyn studied by Kob 
et al. |32| is striking. Furthermore, this change in be- 
havior of X4( t oi)\nve and £ dyn occurs in a temperature 
range where Berthier et al. 29 noticed a change in the 
finite size effects. This change in the finite size effects is 
correlated with an avoided mode-coupling transition. It 
occurs in the temperature range where reasonable mode- 
coupling fits predict the mode-coupling temperature in 
the binary harmonic spheres studied in this, Kob's, and 
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Berthier's work. More work seems to be needed to under- 
stand the complex dynamics around the mode-coupling 
temperature range. 
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As discussed in Sec. |III[ our results for the mode- 
coupling temperature depend on the range of temper- 
atures used for fitting. The fit using the temperature 
range 12 > T > 6 results in T c = 5.1 ± 0.1 which is 
slightly smaller but within error bars of T c = 5.2 of Kob 
et al. [32] . Here we examine some alternative fits that 
can be used for the a-relaxation time and diffusion coef- 
ficient data obtained from our simulations. In particular, 
we show that in addition to the dependence of the mode- 
coupling temperature on the range of temperatures used 
for fitting, there is also a somewhat unexpected depen- 
dence on the quantity being fit. 

In Fig. [15] we show the temperature dependence of the 
a-relaxation time. The dashed lines are fits to a mode- 
coupling like power law r Q = a{T — XjT) -7 . A fit using 
the temperature range 12 > T > 6 results in T C T = 5.2 ± 
0.1 and 7 1 " = 2.7 ± 0.1, but a fit using the range 12 > 
T > 7 results in TJ = 5.8 ± 0.1 and 7 T = 2.3 ± 0.1. Here 
the superscript r denotes the mode-coupling temperature 
obtained from the a-relaxation time fits. We note that 
the former result agrees very well with the result obtained 
by Kob et al. [H [5] by fitting the a-relaxation times 
obtained from the self-intermediate scattering functions 
for small and large particles. 

Also shown in Fig. [15] are two other commonly used 
fits, r a — T Q e k ° T (solid line) and r Q = T v e kv ( T ~ T °* > 
(dotted line). Both fits were performed for T < 8, and 
both fit the data very well for this temperature range. 
The latter fit results in TJ = 2.2 where, again, the super- 
script t denotes the Vogel-Fulcher temperature obtained 
from the a-relaxation time fit. 

In Fig. [16] we show the temperature dependence of 
the diffusion coefficient. The dashed lines are the mode- 
coupling- like fits D = b(T — T®) 1 for the same ranges of 
temperatures as used in the analysis of the a-relaxation 
time. Fit using the range 12 > T > 6 results in 
T ( D = 4.8±0.1 and j D = 2.4±0.1, but a fit for 12 > T > 7 
results in TjP = 5.4±0.1 and 7^ = 2.0±0.1. Here the su- 
perscript D denotes the mode-coupling temperature ob- 
tained from the diffusion coefficient fits. 

We note that unlike what was found in the pioneering 
analysis of a binary Lennard- Jones system by Kob and 



FIG. 15: The a-relaxation time as a function of the inverse 
temperature. The dashed lines are mode-coupling-like fits to 
the function r a = a(T — T^) -7 where the parameters of the 
fits are given in the figure. Two other fits are shown in the 
figure, a solid line which is a fit to r a — roe k ° T for T < 8, 
and a dotted line which is a fit to r a = T„e fcl, ' T-T ° -* with 
To = 2.2. For the latter fits we used T < 8 data. At the 
lowest temperatures the solid and dotted line nearly overlap. 

Andersen [3], we found that the mode-coupling temper- 
atures determined from fits of the a-relaxation time and 
the diffusion coefficient are somewhat different. However, 
we also note that a visual comparison of Fig. |6]with Figs. 
|15| and [TB] shows that in the temperature ranges used for 
fitting the fits imposing the same mode-coupling tem- 
perature for the a-relaxation time and the diffusion co- 
efficient seem as good as the fits allowing for different 
mode-coupling temperatures. 

Also shown in Fig.[l6]are two other commonly used fits, 
D = D e- daT ~ 2 (solid line) and D = B v e- d ^ T - T oY x 
(dotted line). Again, both fits were performed for T < 8, 
and both fit the data very well for this temperature range. 
The latter fit results in = 1.5 which, again, is different 
from the Vogel-Fulcher temperature obtained from the a- 
relaxation time fit. 

To explore the fitting procedure in more depth, we 
also examined the dynamics of the small and large par- 
ticles separately. To this end we calculated F£ (t) and 
<^r 2 ( a )(i)) where the calculations included sums over 
particles of type a, i.e. just the large or just the small 
particles. Then we defined a r° as when F£(t£) = er 1 
and obtained the diffusion coefficients from (5r ( a \t)) — 
6D a t at long times. Then we fitted mode-coupling power 
laws to t° and D a , four sets of data, and imposed that the 
fits have the same mode-coupling temperature T c . First, 
we imposed the restriction that t° for both small and 
large particles have the same exponent j T and that the 
exponent for the diffusion coefficients 7^ are the same, 
but 7 1 " and j D are allowed to be different. For this case, 
using the temperature range 12 > T > 6 we obtained 
7 T = 2.9 ± 0.1 and -y D = 2.3 ± 0.1 with T c = 5.1 ± 0.1. 
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FIG. 16: The inverse of the diffusion coefficient as a func- 
tion of the inverse temperature. The dashed lines are mode- 
coupling-like fits D ~ (T — T,P) 7 with the values of the fit 
parameters given in the figure. Two other fits are shown in 
the figure, a solid line which is a fit to D — Doe~ d ° T and 
a dotted line which is a fit to D = L>„e" a " u ~ J ° ; with 
T° = 1.6. For the latter fits we used T < 8 data. 



Using the temperature range 12 > T > 7 we obtained 
7 r = 2.5±0.1 andj D = 2.1±0.1 with T c = 5.5±0.1. Sec- 
ond, we allowed the exponents to be different for the dif- 
ferent types of particles as well as different for D a and r°. 
For these fits, using the temperature range 12 > T > 6 
we obtained 7J = 2.7 ± 0.1 and = 1.9 ± 0.1 for the 
small particles, = 2.9 ±0.1 and 7 ; D = 2.4 ±0.1 for the 
large particles, with aT c = 5.1 ± 0.1. Using the temper- 
ature range 12 > T > 7 we obtained 7J = 2.3 ± 0.1 and 
7f = 1.7±0.1 for the small particles, H = 2.4±0.1 and 

lP 



yf 

- D = 2.1±0.1 for the large particles, withaT c = 5.7±0.1. 



Ill 



The comparison of the T c values discussed in Sec 
with the T c values presented in this paragraph shows that 
the value of the mode-coupling temperature depends only 
weakly on whether quantities pertaining to all particles 
or quantities pertaining to small and large particles are 
used for fitting. 



To aid in our fits to Si(q; r a ) we calculate the ensemble 
independent dynamic susceptibility 



E 



w n (t)w m (t) 



E W "W 



(11) 

which is the fluctuation in the overlap function ^ w n (t) 
where w n (t) — Q(a — \r n (t) — r n (Q)|) and 0(x) is Heavy- 
side's step function. Since Xi{t) 1S a fluctuation quantity 
and some global fluctuations are suppressed in the NVE 
and NVT simulations presented in this work, one cannot 
calculate Xi{t) directly in an NVE or NVT simulation. 
By using methods developed previously [211 E21 HI] we 
account for these suppressed fluctuations, and this gives 
us an added benefit of examining how dynamic hetero- 
geneities may be related to experimentally determinable 
quantities @1 HI SMS] ■ 

First we consider a transformation between the grand- 
canonical ensemble where the total number of particles 
and the energy are allowed to fluctuate to the canoni- 
cal ensemble where the energy fluctuates, but the total 
number of particles are constant. In what follows we will 



denote 
semble, 
and 



')nvt 



as an average in the grand-canonical en- 
as an average in the canonical ensemble, 



I NVE 



as an average in the micro-canonical ensem- 
ble. Note that there are two types of particles, thus there 
are two chemical potentials [i\ and fi2, but we use a short- 
hand notation where fi represents both \i\ and \ii and N 
represents N\ and N 2 . Thus (•) NVT is an average in an 
ensemble where N\ and N2 are not allowed to fluctuate. 
Using the method of Lebowitz et al. [42; , we obtain 



(N) ^VTXi{t) — Xi{t)\8N + Xi{t)\NVT (12) 



where 



Xi{t)\sN = ((SN,) 2 ) 



fiVT 



NVT 



dNi 

d (En W n) NVT 
8N 2 



-<(<^)Vt 

, & {Hn w n) NVT ® (Sn w n) NVT_ 



dN x 



dN 2 



In conclusion, for the system studied in this paper, we 
found that the mode-coupling temperature seems to be 
somewhat ambiguous. For this reason, in this paper we 
used the notion of the mode-coupling temperature range 
or T c range. On the basis of the fits discussed in Sec. [TTT] 
and in this Appendix, we determined the T c range to be 
5.6-5.1. 



and Xi(t)\NVT is the fluctuation of ^ n w n (t) calculated 
in the canonical ensemble. For the lowest two tempera- 
tures in this work, T = 5, and 5.5 we used Eq. (12) to 
calculate Xi(t)- The temperatures T > 6 are discussed 
further below. 

In practice, we convert the derivative with respect to 
particle numbers to derivative with respect to concen- 
tration c = Ni/N and volume fraction (j> = 4%(Nidf + 
A 2 c?2)/(3U), where d\ = <Ju and d 2 = a 22 . This results 
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in the following expression 

Xi(t)\sN = X%Hi + X4>XcH 2 + Xc#3 

+F {t) 2 H 4 + F {t) X ^H b + F {t) X cH 6 , 

(14) 

where Xx = dF Q (t)/dx. The H n are functions of S nm = 
lim g ^.Q S nm {q) where 



S nm (q) = (iV m JV n )-V 2 



(15) 



are the partial structure factors. The H n terms are given 
by 

Hi = (^) 2 [d 6 1 x 1 S 11 +2d 3 1 d 3 2 ^xTx^S 12 +d 6 2 x 2 S 22 ] 

(16) 



Ho 



n P r^3 



+d 3 ,x 2y /x 1 x 2 S 12 - dlx 1 x 2 S 22 ] (17) 
H 3 = x\x\S\\ - 2x 1 x 2y / x 1 x 2 S 12 + x\x 2 S 22 (18) 
Ha, = x 1 Sn + 2y , x\x 2 Si 2 + x 2 S 22 (19) 



H h 



r^3 



[dlx^u + {d\ + dl)^/xTxiS 12 + dlx 2 S 22 ] 



(20) 



H 6 = 2[xix 2 Sii + (x 2 -xi)t/xJx 2 Si 2 -xix 2 S 22 ], 

(21) 

where p = N/V is the number density, x n — N n /N, and 
we have used that 



WfiVT (SNndNm)^!, = limy/x n x m S nm (q) 



(22) 



In Eqs. (161 through (21 1 the derivative with respect to 



volume fraction is performed at constant concentration 
and the derivative with respect to concentration is per- 
formed at constant volume fraction. One can also per- 
form the calculations varying the number density and 
the concentration instead of the volume fraction and the 
concentration, but we found that this procedure involved 
large cancellations of terms of opposite sign. We also note 
that the Xc term could not be neglected in this work, but 
could be neglected in the analysis of earlier simulations 
of hard sphere systems [22] . 

For T > 6 we performed only NVE simulations, 
thus in Eq. 12 we used XiifilNVT = Xi{t)NVE + 
kBT 2 XT(t) 2 /c v where c v is the constant volume spe- 
cific heat per particle and Xi(t)\NVE is the fluctuation 
of J2 n w n(t) calculated in the NVE simulation. To cal- 
culate all the derivatives we performed at least two sim- 
ulations at x + Sx and two at x — Sx where x is tf>, T, or 
c. The size of Sx depended on the temperature. 

We compared Xii^NVT and Xa^nve + 
kBT 2 XT(t) 2 /c v and found excellent agreement for 



3 50 

X 45 

>40 
Z 

C" 35 
^30 



■ - Z 4 W 

O X 4 (i)I NV t 

— X4» l NVE + X 4 (t)l 8T 
X 4 <t)l NV E 

■ ■ ■ X4WW 



/"\ T=5.5 




FIG. 17: The four-point susceptibility calculated in an NVT 
simulation X4(t)|jvvT (circles) compared to the susceptibility 
calculated in an NVE simulation X4(t)\NVE (dashed line), 
Xi{t)sT (dotted line), and the sum x^)\nve +X4(t)sT (solid 
line). The dashed-dotted line is the full ensemble independent 
susceptibility X-t(t)' 



T = 5.5, see Fig. [17] For comparison we show the 
full ensemble independent susceptibility Xi(^) as the 
dotted-dashed line in Fig. [17] 

We note that good agreement between Xi(t)NVT and 
Xi(t)NVE + kBT 2 xr(t) 2 /c v does depend on the simula- 
tion parameters. Recall that we utilized a Nose-Hoover 
thermostat for the NVT simulations and the LAMMPS 
simulation package [36] which utilizes the equations of 
motion of Shinoda et. al [50]. In the LAMMPS input 
script [371 El] 1 one has to provide a temperature damp- 
ing parameter, Tdamp, which is related to the coupling 
to the heat bath. The Tdamp parameter is roughly equal 
to the time it takes for the temperature to relax. A too 
small Tdamp and the temperature fluctuates wildly, but 
a too large Tdamp results in a large time for the tem- 
perature to equilibrate and a drift in the temperature 
for the very long simulations needed in this work. We 
found that a Tdamp of 2 time units to result in large 
fluctuations, thus a large X<t(t)\NVT- However, notice- 
able temperature drift occurs after about several hun- 
dred million time steps for a Tdamp of 10,000 time units. 
We chose a Tdamp of 1,000 time units to minimize the 
effects of the thermostat on the temperature fluctuations 
and to remove the energy drift. Note that we checked 
that this choice of Tdamp did not influence the results at 
T = 5.5 by comparing simulations of 100t q in the NVE 
and NVT ensembles. We also compared results of NVE 
simulations of 10r Q and NVT simulations of 100r a for 
T = 5 and found statistically the same results. However, 
the statistics are poor for the shorter NVE simulations 
and we needed the longer NVT simulations to obtain 
sufficient statistics for this work. 

We fit Si(q;T a ) using the value of Xi{ T a) calculated 
through the procedure described in this Appendix as 
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FIG. 18: The ensemble independent susceptibility X±(£) f° r 
T — 10, 9, 8, 7, 6, and 5 listed from the smallest peak to the 
highest peak. T hese are the same temperatures as shown for 
S 4 {q;T a ) in Fig. [10 



S4, (q — 0; r a ). This allows us to use slightly smaller sys- 
tems and makes it possible to reach temperatures below 
the T c range. However, to examine what occurs with- 
out utilizing Xi{ T a), we also fitted S^{q; T a ) using only 
finite q data. Thus, we found Xi( T a) purely through the 



extrapolation of Si(q\T a )- These fits generally result in 
a lower value of xa{ t o) and a smaller ^4(r Q ), especially 
when the correlation length ^4(r Q ) and the susceptibility 
Xi( T a) are large enough such that the plateau region of 
the Ornstein-Zernicke function is not visible in the finite 
q data. However, the finite q fits result in values that are 
are always within error of Xi{t) determined from fits us- 
ing also the q = point. Moreover, the correlation length 
^4(t q ) obtained using only finite a data is close to, but 
systematically smaller than £,a{to) determined from fits 
using also the q — point. For most of the temperatures 
the two values of ^4(r a ) are within error bars. 

Shown in Fig. [18] is the full ensemble independent sus- 
ceptibility Xi(t) as a function of time for T = 10, 9, 8, 
7, 6, and 5. At t = Xi{t) coincides with the small 
wave-vector limit of the total static structure factor, 



Xi {t = 0) 



lim N~ 



-jq-[r„(0)-r ro (0)] 



lim S(q). 

\m,Ti 1 

(23) 

Xi(t) begins to grow during the /3- relaxation and it peaks 
around r a . For longer times Xi(t) decays to zero. The 
growth towards the peak can be well described by a power 
law Xi(t) = X + xot c where c depends on temperature 
and decreases with decreasing temperature in the tem- 
perature range studied. 
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